/***********************************************************************
Compile       : gcc -O3 -o sequential sequential.c ../design/design.c ../design/utils.c -I../design
Usage         : ./sequential mode mode1 mode2 irflag N L inputfile 
                (./sequential 0 0 0 1 11 12 Budding.txt)
mode          : 0 normal ; 1 cycled ; 2 triggled cycled
mode1         : 0 disallow self-activation ; 1 allow self_activation
mode2         : 0 no Restrictions ; 1 Restriction I ; 2 Restriction II; 3 Restriction I and II
irflag        : 1 irreducible solution only
N             : number of the nodes
L             : number of the steps        
Output        : activity min_link designability1 designability2
activity      : number of 1 in the process
min_link      : minimum links of the qualified network
designability1: the minimum designability for one node
designability2: total designability
***********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include "design.h"

int counter1, counter2, counter3, counter4;

void print_help ()
{
	printf ("Usage         : ./sequential mode mode1 mode2 irflag N L inputfile\n");
	printf ("                (./sequential 0 0 0 1 11 12 Budding.txt)\n");
	printf ("mode          : 0 normal ; 1 cycled ; 2 triggled cycled\n");
	printf ("mode1         : 0 disallow self-activation ; 1 allow self_activation\n");
	printf ("mode2         : 0 no Restrictions ; 1 Restriction I ; 2 Restriction II; 3 Restriction I and II\n");
	printf ("irflag        : 1 irreducible solution only\n");
	printf ("N             : number of the nodes\n");
	printf ("L             : number of the steps\n\n");
	printf ("Output        : activity min_link designability1 designability2\n");
	printf ("activity      : number of 1 in the process\n");
	printf ("min_link      : minimum links of the qualified network\n");
	printf ("designability1: the minimum designability for one node\n");
	printf ("designability2: total designability\n");
}

void init_design (struct designability *design, int argc, char *argv[])
{
	int N, L, i;

	design->mode = atoi (argv[1]);
	design->mode1 = atoi (argv[2]);
	design->mode2 = atoi (argv[3]);
	design->irflag = atoi (argv[4]);
	design->bdflag = 0;
	design->dgflag = 0;
	design->pnflag = 1;
	design->svflag = 1;
	design->N = atoi (argv[5]);
	design->L = atoi (argv[6]);

	design->N += design->mode1;
	N = design->N;
	L = design->L;

	/* the degree of irreducible solution is zero */
	if (design->irflag)
		design->dgflag = 0;

	design->pathway = calloc (sizeof (int), N * L);
	design->NS2 = calloc (sizeof (int), N);
	design->min_link = calloc (sizeof (int), N);
	design->links_dist = calloc (sizeof (int), N * (N + 1));
	design->Links_dist = calloc (sizeof (double), N * N + 1);
	design->degree_dist = calloc (sizeof (int), N * (2 * N + 1));
	design->Degree_dist = calloc (sizeof (double), 2 * N * N + 1);
	design->solution = calloc (sizeof (struct solution), N);

	for (i = 0; i < N; i++)
	{
		init_solution (&(design->solution[i]));
	}
}

void free_design (struct designability *design)
{
	int i;

	free (design->pathway);
	free (design->NS2);
	free (design->min_link);
	free (design->links_dist);
	free (design->Links_dist);
	free (design->degree_dist);
	free (design->Degree_dist);
	free (design->solution);

	for (i = 0; i < design->N; i++)
	{
		free_solution (&(design->solution[i]));
	}
}

void find_connect (int N, int network[N][N], int flag[N], int i)
{
	int j;

	for (j = 0; j < N; j++)
	{
		if (flag[j] != -1)
			continue;
		if (network[j][i] == 1 || network[i][j] == 1)
		{
			flag[j] = flag[i];
			find_connect (N, network, flag, j);
		}
	}
}

/* Return Value : number of subgraphs */
int calc_subgraphs (int N, int network[N][N], int flag[N])
{
	int i, num;

	for (i = 0; i < N; i++)
		flag[i] = -1;

	num = 0;
	for (i = 0; i < N; i++)
	{
		if (flag[i] == -1)
		{
			flag[i] = num;
			num ++;
			find_connect (N, network, flag, i);
		}
	}

	return num;
}

/* Return Value: 1 sequential 0 non-sequential */
int check_sequential1 (int N, int network[N][N])
{
	int i, j, p, flag;

	for (i = 0; i < N; i++)
	{
		flag = 0;
		for (j = 0; j < N; j++)
		{
			if (i == j)
				continue;

			if (network[i][j] == 1)
			{
				flag = 1;
				break;
			}
		}
		if (flag == 0)
			break;
	}

	if (flag == 1)
		return 0;
	else if (N == 2)
		return 1;
	else
	{
		int x, y, network2[N - 1][N - 1];

		p = i;

		for (i = 0; i < N; i++)
		{
			if (i == p)
				continue;
			else if (i > p)
				x = i - 1;
			else
				x = i;
			for (j = 0; j < N; j++)
			{
				if (j == p)
					continue;
				else if (j > p)
					y = j - 1;
				else
					y = j;

				network2[x][y] = network[i][j];
			}
		}
		return check_sequential1 (N - 1, network2);
	}
}

void setup_level (int N, int network[N][N], int level[N])
{
	int i, j, flag, curr_level, n;

	/*
	printf ("Parsing network:\n");
	for (i = 0; i < N; i++)
	{
		for (j = 0; j < N; j++)
		{
			if (network[i][j] == 1)
				printf ("1 ");
			else
				printf ("0 ");
		}
		printf ("\n");
	}
	*/

	curr_level = n = 0;
	for (i = 0; i < N; i++)
		level[i] = -1;

	while (1)
	{
		for (i = 0; i < N; i++)
		{
			if (level[i] != -1)
				continue;
        
			flag = 0;
			for (j = 0; j < N; j++)
			{
				if (i == j)
					continue;
        
				if (level[j] >=0 && level[j] < curr_level)
					continue;
        
				if (network[i][j] == 1)
				{
					flag = 1;
					break;
				}
			}
			if (flag == 0)
			{
				level[i] = curr_level;
				n ++;
				//printf ("Node %2d Level %2d N %2d\n", i, curr_level, n);
			}
		}
		if (n == N)
			break;
		curr_level ++;
	}
}

/* Return Value: 1 sequential (more restrictive definition) */
int check_sequential2 (int N, int network[N][N])
{
	int i, j, level[N], flag;
	
	setup_level (N, network, level);

	flag = 1;
	for (i = 0; i < N; i++)
	{
		for (j = 0; j < N; j++)
		{
			if (network[i][j] == 1)
			{
				if (level[i] != level[j] + 1)
				{
					flag = 0;
					break;
				}
			}
		}
		if (flag == 0)
			break;
	}

	return flag;
}

/* Return Value: 1 sequential (sequential1[]: weak; sequential2[]: strong)*/
void check_sequential (int N, int network[N][N], int subgraph[N], int M, int sequential1[M], int sequential2[M])
{
	int i, j, k, n;

	for (i = 0; i < M; i++)
	{
		n = 0;
		for (j = 0; j < N; j++)
		{
			if (subgraph[j] == i)
				n ++;
		}

		if (n == 1)
		{
			sequential1[i] = 1;
			sequential2[i] = 1;
			continue;
		}

		int network2[n][n], idx[n];

		n = 0;
		for (j = 0; j < N; j++)
		{
			if (subgraph[j] == i)
			{
				idx[n] = j;
				n ++;
			}
		}

		for (j = 0; j < n; j++)
		{
			for (k = 0; k < n; k++)
			{
				network2[j][k] = network[idx[j]][idx[k]];
			}
		}

		sequential1[i] = check_sequential1 (n, network2);

		if (sequential1[i] == 1)
			sequential2[i] = check_sequential2 (n, network2);
		else
			sequential2[i] = 0;
	}
}

void print_complete_solution (struct designability *design, int N, int select[N], int node)
{
	int i, j, k;

	if (node == N)
	{
		counter1 ++;
		/*
		printf ("Solution %d\n", counter1);
		for (i = design->mode1; i < N; i++)
		{
			if (i == design->mode1 && design->mode == 2)
				continue;

			j = select[i];
			for (k = design->mode1; k < N; k++)
			{
				if (k == design->mode1 && design->mode == 2)
					continue;

				if (design->mode1 == 1 && k == i && design->solution[i].edge[j * N] == 1)
					printf (" g");
				else if (design->solution[i].edge[j * N + k] == 0)
					printf (" n");
				else if (design->solution[i].edge[j * N + k] == 1)
					printf (" g");
				else if (design->solution[i].edge[j * N + k] == -1)
					printf (" r");
			}
			printf ("\n");
		}
		*/

		int M;

		M = design->mode1;

		if (design->mode == 2)
			M ++;

		int network[N - M][N - M];

		for (i = design->mode1; i < N; i++)
		{
			if (i == design->mode1 && design->mode == 2)
				continue;

			j = select[i];
			for (k = design->mode1; k < N; k++)
			{
				if (k == design->mode1 && design->mode == 2)
					continue;

				if (design->mode1 == 1 && k == i && design->solution[i].edge[j * N] == 1)
					network[i - M][k - M] = 1;
				else if (design->solution[i].edge[j * N + k] == 0)
					network[i - M][k - M] = 0;
				else if (design->solution[i].edge[j * N + k] == 1)
					network[i - M][k - M] = 1;
				else if (design->solution[i].edge[j * N + k] == -1)
					network[i - M][k - M] = -1;
			}
		}

		int subgraph[N - M], n;

		n = calc_subgraphs (N - M, network, subgraph);
		if (n == 1)
		{
			counter3 ++;
		}

		int sequential1[n], sequential2[n];

		check_sequential (N - M, network, subgraph, n, sequential1, sequential2);

		int flag1, flag2;

		flag1 = flag2 = 0;
		for (i = 0; i < n; i++)
		{
			flag1 += sequential1[i];
			flag2 += sequential2[i];
		}

		if (flag1 == n)
			counter2 ++;
		if (flag2 == n)
			counter4 ++;
	}
	else
	{
		for (i = 0; i < design->solution[node].curr; i++)
		{
			select[node] = i;
			print_complete_solution (design, N, select, node + 1);
		}
	}
}

void print_solution (struct designability *design)
{
	int i, j, k, N;

	N = design->N;
	if (design->svflag && design->pnflag)
	{
		for (i = design->mode1; i < N; i++)
		{
			printf ("Total Number of Soultion for Node %2d : %d\n", i - design->mode1, design->NS2[i]);
			for (j = 0; j < design->solution[i].curr; j++)
			{
				printf ("Node %2d :", i - design->mode1);
				for (k = design->mode1; k < N; k++)
				{
					if (design->mode1 == 1 && k == i && design->solution[i].edge[j * N] == 1)
						printf (" g");
					else if (design->solution[i].edge[j * N + k] == 0)
						printf (" n");
					else if (design->solution[i].edge[j * N + k] == 1)
						printf (" g");
					else if (design->solution[i].edge[j * N + k] == -1)
						printf (" r");
				}
				printf ("\n");
			}
		}
	}
}

void print_pathway (struct designability *design)
{
	int i, j;

	printf ("c %3d %3d %5d %lg\n", design->activity, design->min_links, design->NS2_min, design->NS);

	for (i = 0; i < design->L; i++)
	{
		for (j = design->mode1; j < design->N; j++)
		{
			printf ("%d ", design->pathway[i * design->N + j]);
		}
		printf ("\n");
	}
}

int main (int argc, char *argv[])
{
	struct designability design;

	if (argc == 1)
	{
		print_help ();
		exit (0);
	}

	init_design (&design, argc, argv);

	FILE *fp;
	char buf[256];
	int i, j, N, L;

	N = design.N;
	L = design.L;

	fp = fopen (argv[7], "r");

	while (fgets (buf, 256, fp))
	{
		design.activity = 0;
		for (i = 0; i < L; i++)
		{
			for (j = 0; j < N; j++)
			{
				if (design.mode1 == 1 && j == 0)
				{
					design.pathway[i * N + j] = 1;
					continue;
				}
				fscanf (fp, "%d", &(design.pathway[i * N + j]));
				if (i <= L - 3 || (i == L - 2 && design.mode != 2))
					design.activity += design.pathway[i * N + j];
			}
			fgets (buf, 256, fp);
		}

		if (design.mode == 2)
			design.activity --;

		/* Get the number of solution using enumerating algorithm */
		if (calc_design (&design) == 0)
			continue;

		print_pathway (&design);
		//print_solution (&design);

		if (design.svflag && design.pnflag)
		{
			int select[N];

			counter1 = counter2 = counter3 = counter4 = 0;

			if (design.mode == 2)
				print_complete_solution (&design, N, select, design.mode1 + 1);
			else
				print_complete_solution (&design, N, select, design.mode1);

			printf ("Single Connected Graph : %4d / %4d\n", counter3, counter1);
			printf ("Check Sequentialilty : %4d %4d / %4d\n", counter4, counter2, counter1);
		}
	}

	fclose (fp);

	free_design (&design);

	exit (0);
}
